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ABSTRACT 

We obtain the non-linear relation between cosmological density and velocity per- 
turbations by examining their joint dynamics in a two dimensional density-velocity 
divergence phase space. We restrict to spatially flat cosmologies consisting of pres- 
sureless matter and non-clustering dark energy characterised by a constant equation 
of state w. Using the spherical top-hat model, we derive the coupled equations that 
govern the joint evolution of density and velocity perturbations and examine the flow 
generated by this system. In general, the initial density and velocity are independent, 
but requiring that the perturbations vanish at the big bang time sets a relation be- 
tween the two. This traces out a curve in the instantaneous phase space, which we 
call the 'Zeldovich curve'. We show that this curve acts like an attracting solution 
for the phase space dynamics and is the desired non-linear extension of the density- 
velocity divergence relation. We obtain a fitting formula for the curve as a function 
of il m and w and find that, as in the linear regime, the explicit dependence on the 
dark energy parameters stays weak even in the non-linear regime. Although the result 
itself is somewhat expected, the new feature of this work is the interpretation of the 
relation in the phase space picture and the generality of the method. Finally, as an 
observational implication, we examine the evolution of galaxy cluster profiles using 
the spherical infall model for different values of w. We demonstrate that using only 
the density or the velocity information to constrain w is subject to degeneracies in 
other parameters such as erg but plotting observations onto the joint density-velocity 
phase space can help lift this degeneracy. 

Key words: cosmology: theory - cosmology: dark energy - cosmology: large-scale 
structure of Universe - galaxies: clusters: general 



INTRODUCTION 



The large scale distribution of matter into overdense and underdense regions and their associated peculiar velocities are the 
two main observables that indicate that our universe is far from homogeneous. They are related via the continuity equation, 
the Euler equation and the law of gravitation. In the linear regime, the relation between the peculiar velocity 'v' and the 
fractional overdensity '<5' is V • v = —fHS, where H is the Hubble parameter and / is the linear growth factor. / mainly 
depends on the matter density parameter Q m , and is usually expressed as QZi, where 7 is the growth index. Given 7, a 
comparison of the data from redshift and peculiar velocity surveys constrains fi m , or some combination of Q m and the galaxy 
bias parameter. Early studies in this field mostly assumed a pure matter universe and used 7 ~ 0.6 ( Peebles||1976[ ), to either 
get bias independent measures of mass from velocity fields or to constrain Sl m (see review articles by Dekel|l9*9*4 



and Strauss 



fc Willick||1995 1. The estimate by Peebles, given more than thirty five years ago, turned out to be a good approximation for 
a large range of models. The dependence of the linear growth rate on the cosmological constant (for e.g., |Martel|1991 Lahav 
et al.|1991 | or on the dark energy equation of state w | |Wang fc Steinhardt|1998]|Linder|2005[ ) was shown to be relatively weak. 
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However, more recently, it has been pointed out that models with similar expansion histories but with different gravitational 
dynamics can be distinguished by their growth indices; 7 is 0.55 for the ACDM model, but 0.67 for certain modified theories 
of gravity (Lue, Scoccimarro, & Starkman 20041. Accordingly, modern observational efforts are focussed on measuring the 
linear growth rate with the added aim of constraining 7 (for e.g., Guzzo et al.|2008 Majerotto et al.|2012 |. 

The method of velocity reconstruction from redshift surveys using linear theory breaks down when the perturbations 
are of order unity. Quasi-linear effects need to be included even to get an accurate determination of the linear growth rate 
(Nusser, Branchini, & Davis 2012). Various analytic quasi-linear and non-linear extensions using second order perturbation 



theory (Bernardeau 



1992 



hereafter B92), the Zeldovich approximation (Nusser et al. [19911 and higher order Lagrangian 



perturbation theory (Chodorowski & Lokas 1997; Chodorowski et al.|1998 |Kudlicki et al.|2000[ ) have been proposed. Accurate 
numerical results can be obtained using N-body simulations, however, N-body codes give mass-weighted instead of volume- 
weighted estimates ( Dekel, Bertschinger, fc Faber|1990 l making it difficult to compare them with analytical answers. Refined 
volume averaged estimates were given using uniform-grid codes ( Kudlicki et al.p OOO ; Ciccielag ct al. 2003 1 and the method 
of tessellations ( Bernardeau fc van de Weygaert|1996 Bernardeau et al.|1999 k In general, these simulations are slow and are 
applicable to a limited range of models. It is usually assumed that the weak dark energy dependence of the linear relation 
extends to the non-linear regime and often the results are quoted in terms of the scaled velocity divergence V • v//(f2 m ), 
absorbing the cosmology dependence into the linear growth rate. While this assumption has been tested for different values 
of the cosmological constant A (for e.g., Lahav et al.||199~ I , the explicit dependence on the dark energy equation of state w 
has not yet been derived. 

The spherical top-hat system is an alternate way to model evolution in the non-linear regime. The model is simple to 
solve; equations of motion reduce to ordinary second order differential equations for the evolution of the scale factor. The main 
drawback is that it does not take into account interaction between scales; the price paid for computational ease. Nevertheless, 
it has been successful in predicting non-linear growth until virialization ( Engineer et al.|20 00 Shaw fc Mota|2008 k Recently 



Bilicki & Chodorowski (20081, hereafter BC08, obtained the non-linear density- velocity relation from the spherical top-hat. 



They utilised the top-hat's known exact analytic solutions and hence were restricted to pure matter cosmologies. We adopt a 
different approach, one that can be generalised to a range of background cosmologies. We numerically investigate how generic 
density and velocity perturbations evolve in a two dimensional density-velocity divergence phase space. We identify a special 
curve which is obtained by imposing the condition that the perturbations vanish at the big bang time and show that it acts 
like an attracting solution for the dynamics of the system. We refer to this as the 'Zeldovich curve' and demonstrate that 
it is the desired density-velocity relation in the non-linear regime. This approach was first put forward in a recent paper 
( Nadkarni-Ghosh fc Chernoff|20lT hereafter NC), but the analysis there was restricted to a pure matter, Q. m = 1 cosmology. 
Here, we extend the same idea to flat cosmological models with pressureless matter and non-clustering dark energy described 
by a constant equation of state w. 

The paper is organised as follows. [j2]gives the equations governing the spherical top hat, introduces the 'Zeldovich curve' 
and demonstrates its importance in the evolution of perturbations in the joint density-velocity phase space. !^4] defines the 
non-linear growth factor and discusses its dependence on the dark energy equation of state. Sj3] gives a fit to the Zeldovich 
curve and compares it to others in the literature. Sj5]discusses the evolution of galaxy cluster profiles using the spherical infall 
model and examines this curve in the context of observationally relevant quantities. Sj6] gives the summary and conclusion. 



2 DYNAMICS OF THE SPHERICAL TOP-HAT 
2.1 Physical set-up and equations 

The compensated spherical top-hat perturbation consists of a uniform density sphere surrounded by a concentric spherical 
compensating region; homogenous and isotropic background extends beyond the outer edge of this region. The compensating 
shell could be empty or contain mass depending upon whether the inner sphere is overdense or underdense with respect 
to the background. The background cosmology is flat and comprises of pressureless matter and dark energy described by 
a constant equation of state parameter w, not necessarily —1. Only matter perturbations are considered. Dark energy is 
assumed to be spatially uniform and does not interact with the matter fields. The total matter density and Hubble parameter 
are denoted as p m (pm) and H(H) for the background (perturbation). The dark energy density is the same for the background 
and perturbation and is denoted as p<j,. 

The evolution of the background is completely determined by specifying a^, Hi p m ,i and p<£,i at some arbitrary initial 
time ti. The governing equation is 

a m 



a 2 



^ + (1 + 3^,^) 



3(1+™) 



(1) 



where the dots represent derivatives with respect to time t and Q m (<t,),i are the initial matter (dark energy) density parameters. 
The initial conditions for eq. |l]) are a(ti) = a,i and a(ti) — a,iHi. 
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Figure 1. Variation of f2 m (a) for a flat cosmological model with dark matter and dark energy described by a constant equation of state 
w. For all models, Q m ,o = 0.29 and Cl^fi = 0.71. The values coincide at very early times and at the present epoch, but the evolution 
through intermediate epochs is different for different values of w. 



The perturbation at the initial time is described by two quantities: the fractional overdensity Si and the fractional Hubble 
parameter S v ,i defined as 



Si = 



pm . i 
pm,i 

Hi _ 
Hi 



(2) 



(3) 



In the case of a pure matter cosmology, Birkoff's theorem guarantees that the inner sphere can be described as a separate 
independent universe, whose scale factor obeys an equation analogous to eq. but with a different value for the matter 
density parameter. In the presence of a dark energy term the modification is not always obvious as has been discussed by 
various authors ( |Pace, Waizmann, fc Bartelmann|2010| [Wintergerst fc Pettorino|2010[ ). We follow the approach advocated in 
these papers of starting with the exact non-linear hydrodynamic equations for the density and velocity instead of an equation 



for the perturbation scale factor (for e.g., Lima, Zanchin, & Brandenberger 1997 Abramo et al. 20071. By recasting the 
Eulerian system in Lagrangian coordinates, one can show that for scales much smaller than the horizon scale and in the 
absence of matter-dark energy interactions, the spherical top-hat can be described by a scale factor b(t) which obeys 



Hi 
2 



n m ,ia?(l + &) 
b 3 



+ (l + 3w)ft^i (^-) 



3(l+m) 



(4) 



with initial conditions b(ti) = ai and b(ti) = cti ( 1 + 5 Vt i) (see Appendix |a| for details). Hi, Q, m .i, Q<t>,i in eqs. |l]) and Q are 
specified at some arbitrary initial time ti, which is different from the time to denoting 'today' when a(to) = 1. The Hubble 
parameter H and the density parameters at any arbitrary epoch are related to those at to through 



H 2 (a) 



H m (a) = 



Hi 



+ 



ra.O 



a 3(l+w) 
3 



f^0(a) 



.on 



-3 + fV oa -3(l+») ' 

1 „-3(!+™) 
'</>,0 a 



-3 



-3(1+uj) ' 



(5) 
(6) 
(7) 



The last two equalities follow from the definition of fl, the relation between H and the critical density p c and the scalings 



-3(1+™) 



The evolution of Q. m (a) is plotted in figure [l] for four different values of the equation of state 



Flatness implies Q,p(a) = 1 — f2 m (a) at all epochs. f2 m .o and Q$,o are chosen to be 0.29 and 0.71 respectively. All models are 
matter dominated at early epochs (Q m (a « 1) = 1) and have the same matter density fi m ,o = 0.29 at late epochs. But at 
intermediate epochs, f2 m (a) is different for the four models. Dark energy dominates earlier for larger values of w. For a fiat 
cosmology consisting of pure matter, fi m = 1 and = for all epochs. 

Once the solutions to eqs. |l]) and Q are known, the perturbation parameters at any later time can be computed 

(1 + S t )a 3 



6(t) 



6 v (t) = 



1_& 
Hb 



(8) 
(9) 
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Figure 2. Zeldovich curves for different cosmologies at different epochs. Four cosmological models with w = —1/2, —2/3, — 1, —1.5 are 
shown (bottom to top, colour code same as that in figure [l}. The black curve is for the Q m = 1 case; it is the same at all epochs. In 
general, for a fixed value of to, the curves evolve with time. At early times, the various curves coincide with the Q m = 1 and almost 
coincide with each other at the present epoch, deviating at intermediate epochs. This indicates that they are mainly dependent on Q m 
and the explicit dependence on w is very weak. 



The perturbation variable 8 V is related to the usual peculiar velocity divergence as 

Sv = Jh v ' v - (10) 

In the rest of the paper we refer to 8 as simply 'density' and peculiar velocity as simply 'velocity'. 



2.2 The Zeldovich curve 

Consider the two-dimensional phase space whose abscissa and ordinate are 8 and S v respectively. In general, the initial 8 and 
8 V are independent parameters, however, for cosmologically relevant initial conditions, it is usually assumed that there are 
no perturbations at the initial time. This sets a specific relationship between the two. In other words, given an initial 8, S v 
is chosen such that the background and perturbation have the same age. For any starting epoch, this relation traces out a 
curve in the 8 — S v phase space. The 'equal age' condition is imposed in linear theory by ignoring the decaying modes in the 
solution for 8. This assumption also forms the basis of the Zeldovich approximation ( ZerDovich]|1970 1 and hence we refer to 
this special curve as the 'Zeldovich curve'. 

In NC this curve was examined for pure matter cosmologies with fi m = 1. In this special case, from eqs. 

and Q, one 

can explicitly calculate the ages of the background and perturbation and the 'equal age' condition becomes 

^ d JL = 2 (H) 

=0 [(1 + <%- 1 + (1 + <5 u ) 2 -(1 + <5)] 1/2 3 

This relation does not involve any time dependent quantities and the Zeldovich curve has the same form at all epochs. 
However, in the presence of generic dark energy terms, this relation has a time dependence. This can be seen by scaling a and 
b in eqs. and Q by a; and replacing time by tHi. The equations and initial conditions depend only on the parameters 
8i, S Vl i, Q m ,i and w (flatness removes fi^j dependence). The 'equal age' criterion remains unchanged by this scaling and the 
resulting Si — S V: i relation depends only on the parameters Q m ,i and w. Since the value of Q, mi i changes with the starting 
epoch, the relation will be different at different epochs. For the kind of dark energy models considered here, analytic solutions 
for eqs. Q and Q are somewhat tedious, if at all possible ( |Lee fc Ng|2010| ), and we solve the equations numerically. Appendix 
[B] gives the details of the calculation. 

Figure [5] shows the Zeldovich curves for different cosmological models. Brown (long dashed), red (dotted), blue (short 
dashed) and green (dot dashed) lines correspond to fl m ,o = 0.29, £lq>-0 = 0.71 and w = —1.5, — 1, —3/4, —2/3 respectively. The 



black curve, given by eq. ( 11 1, corresponds to the S7 m = 1, Q,^ = model. Six different epochs are considered. At recombination 
(a ~ 0.001), the curves for different w overlap with each other and with the Q m = 1 case. At a = 1 they again overlap, but 
differ from the fi m = 1 case. At intermediate epochs, they differ from each other. Comparison with figure [l] suggests that the 
Zeldovich curve mainly depends on Q. m and the explicit dependence on w is very weak. 

In the next section we will show that these curves play a special role in the dynamics of the perturbations in the 8 — S v 
phase space. 
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Figure 3. Phase space flow for two different cosmological models. The upper panel is the Q m = f model; the Zeldovich curve stays 
the same at all times and is an attractor of the flow. The lower panel is the ACDM model with Q m fi = 0.29, Q<b,a = 0.71, w = — 1. 
The Zeldovich curve is different at different times and so is the instantaneous phase portrait, a signature of the non-autonomous nature 
of the system. The coloured ring of initial perturbations (coloured dots) starts near the origin at a = 0.1. As the evolution proceeds 
the perturbations asymptotically approach the Zeldovich curve (numerical values in text). The curve acts like an attractor and is the 
required non-linear extension of the density-velocity divergence relation. 



2.3 Dynamics in the 5 — 8 V phase space 

The definitions in eqs. |8f and Q combined with eqs. Q and Q give equations that govern the phase space evolution of S 
and S v : 

5 = -3H6 v (l + 5) (12) 
L = -~ [n m (a)5-^{O m (a) + (l + 3w)fi^(a)-2} + 2^] . (13) 

In general, because of the time variation of fl m (a), the system defined by eqs. ( |12[ ) and ( |13[ ) is non-autonomous for dark energy 
models and techniques of linear stability analysis that are usually performed on autonomous systems are not applicable. 
However, it is instructive to first consider the Q. m = 1 case for which the temporal dependence drops and the system becomes 
autonomous. Figure [3] shows the phase portrait for the two qualitatively different cases: Q m = 1 (upper panel) and ACDM 
(lower panel) at three different epochs a — 0.1, 0.6, 1. The blue lines with arrows are the streamlines of the flow. A streamline 
is drawn such that the tangent at any point gives the direction of the phase space velocity vector (8, S v ) at that point. It 
denotes the instantaneous direction of evolution of the system. For an autonomous system, the velocity vector is independent 
of time and the pattern of streamlines is the same at all epochs. However, for a non-autonomous system, the phase space 
velocity has a time dependence and hence the pattern of streamlines changes with a. 

For the tt m = 1 case, the system has three fixed points, shown by the black dots, at (0,0), (—1, —1) and (—1,0.5). The 
saddle point at the origin (0,0) corresponds to a unperturbed cosmology, the unstable node at (—1,-1) corresponds to a 
vacuum static model and the attracting point at (—1,0.5) corresponds to a expanding void model. The solid blue curve is 
the Zeldovich curve. The multi-coloured ring of points at a — 0.1 corresponds to initial conditions with different values of 8 
and 8 V , all at a distance of 0.11 from the centre. Some are very close to the Zeldovich curve at the start, but others are off 
the curve. However, as the flow proceeds, all the points evolve to lie close to the Zeldovich curve at a — 0.6 and a = 1. The 
relative deviations from the curve are estimated as |<5 e „ i. (8 V )/Szei. {S v ) — 1|, where at a given value of 5 V , S evo i. is the value 
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evolved from initial conditions at a = 0.1 and Sz e i. is the value given by the definition of the Zeldovich curve. The closeness 
to the curve is characterised by the maximum relative deviation over the points considered for the evolution. Thirteen points 
were considered and the maximum relative deviation was 14% at a = 0.6 and 3% at a = 1. The deviation decreases with 
time indicating that the trajectories asymptotically approach the curve. A point that starts along the curve at a — 0.1 was 
found to remain along the curve with a maximum relative deviation of 10 -8 over the entire range of evolution from a — 0.1 
to a = 1. Thus, the Zeldovich curve forms an invariant set and acts like an attracting solution for the system. 

For the ACDM case, at a — 0.1, fi m w 1 and Q$ w 0. Hence the flow pattern resembles that of the Q. m = 1 cosmology and 
the Zeldovich curves of the two cases are coincident. But as was discussed in the previous section, the Zeldovich curve evolves 
with time. It is possible to discuss the stability of this system in more formal terms using tools from dynamical systems theory, 
but for the purposes of this paper, we will simply illustrate the 'attracting' behaviour of the Zeldovich curve. We consider 
the same set of initial conditions, represented by the multi-coloured ring a — 0.1, and evolve them using eqs. ( 12 I and (131. 
At a — 0.6 and at a — 1 the set has evolved to lie close the corresponding Zeldovich curves at that epoch. The maximum 
relative deviations from the corresponding Zeldovich curve were defined similar to the earlier case, but at fixed 8 instead of 
5 V . The maximum relative deviation was 17% at a = 0.6 and 6% at a — 1, again indicating that the Zeldovich curve acts like 
an attractor. A point that starts along the curve at a = 0.1 stays along it with a maximum relative deviation of 10 -6 over 
the entire range of evolution. Hence, even in the ACDM case, initial conditions that satisfy the Zeldovich relation at the start 
continue to maintain it and those that do not, evolve such that they establish it. 

Thus, the Zeldovich curve has a special significance in the dynamics of the perturbations as they evolve in phase space. 
It gives the long term behaviour of the density and velocity perturbations and is the desired exact non-linear extension of 
the linear density-velocity divergence relation. Fitting forms for this relation are given in the next section. The fits found are 
accurate to 3%, whereas the deviations measured here were greater than this accuracy and one may argue that the fitting 
forms cannot be useful approximation of the dynamics. However, the initial values of a — 0.1 and the ring of radius 0.11 were 
chosen merely for the sake of convenient plotting and to demonstrate the asymptotic behaviour. Real cosmological initial 
conditions start evolving at recombination (a = 0.001) and sense the presence of the Zeldovich attractor for a longer period 
of time. We evolved a similar ring of initial points with an amplitude of 1.1 x 10 -3 from a = 0.001 to a — 1. The maximum 
relative deviation over the entire range was within 0.01% indicating that the fits in the following section are indeed a good 
approximation to the dynamics at the percent level. 



3 FITTING FUNCTIONS AND COMPARISON WITH EXISTING FORMS 

Previous analytic estimates of the non-linear density-velocity divergence relation have been obtained by B92 using higher 
order perturbation theory and by BC08 using the exact solutions of the spherical top-hat model. Both these estimates were 
for pure matter cosmologies. This section compares these functional forms with the Zeldovich curves obtained from the phase 
space dynamics. We find that a combination of these forms with slight modifications provides a good fit for the curves. 

As was discussed in §2.2| for the flat models considered in this paper, we expect the curves to be functions only of f2 m 
and w. For a given w, Q m varies with a and this corresponds to a family of curves for each value of w. To estimate the 
dependence on the two parameters, Zeldovich curves were obtained for seven different values of w in the range —1/2 to —3/2 
and at eleven epochs in the range a = 0.001 to a — 1. The 8 range for each curve was restricted to —1 ^ 8 ^ 10. For the fit 
in the quasi-linear regime, we chose a generalised form of the B92 formula 

S v = A(Q m ,w)(l-(l + 8) B{nm ' w) ). (14) 

The numerical 8 = — 1 limit of the Zeldovich curve and the slope of the curve in the regime \8\ < 0.01 fixed the two unknowns 
A(Q m ,w) and B(Q m ,w). 

A(Q m ,w) = (15) 

B(n m ,w) = |n^w, (16) 

where 

7i(w) = 0.56(-w)-°' 08 (17) 
72 (u>) = -O.Ol(-w)" 1 ' 18 . (18) 

This fit gave a maximum relative error of 3% in the range —1 ^ 5 ^ 1. It gave a 11% error in a higher range of 8 values 
1 ^ 8 ^ 10. In this range, a modification of the formula of BC08 for positive <5 provides a better fit with 3% maximum relative 
error 

S v = + 8) 1/6 - (1 + <5) 1/2 ), (19) 
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Figure 4. Comparison of the various approximations of the density- velocity relationship in literature for w = — 1 at a = 1. The left 
(right) panels show the relationship for overdensities (voids). The combined 3% fit is a modified version of B92 over the range — 1 <5 SC 1 
and of BC08 over the range 1 < S ^ 10 and does marginally better than the unmodified formulas which were derived for pure matter 
models. 



with the same 71 (w) and 72 (w) defined in eqs. (171 and (18 1. In the linear limit, this fit reduces to 



S„ — — 



n 



71M+72M 



(20) 



For a ACDM cosmology, the growth index 



1998 Linder 20051. This fit also modifies the 



This limit is often important when fits to 



which gives a linear growth index 7 = 71 + 72 = 0.56(— w) 08 — 0.01(— w) 1 
reduces to 7 = 0.55, in agreement with other estimates (Wang & Steinhardt 
void limit (8 = —1) from 8 V — — — g 4 - for the B92 formula to 5 V = — Q.m 56t -~ w ^ 
simulation results are compared to analytical estimates (for e.g., Bernardeau et aL"||1997 Kudlicki et al. [2000 1. 

Figure [4] compares the unmodified forms of B92 and BC08 with the fit presented in this paper for one particular case 
w = — 1 at a = 1. The left panel compares positive overdensities, the right panel compares voids. The original B92 (BC08) fits 
are about 11% (10%) accurate over the entire range —1 ^ 5 ^ 10, implying that the modified version is marginally better. We 
have assumed spherical symmetry to derive the non-linear relation. As was commented in BC08, better agreement with the 
spherical system does not guarantee better agreement with the real system. Numerical simulations may be required to give 
more refined estimates of this relation, nevertheless, Zeldovich curve obtained for the spherical system can provide a starting 
point to analyse the results. 



4 NON-LINEAR GROWTH RATE 

The Zeldovich curve gives the instantaneous location of the perturbations in phase space. It does not directly give information 
about the rate at which perturbations evolve along the curve. This is encoded in the growth rate, usually denned as dlnS/dlna. 
From eq. (1 1 2h . we see that the growth rate for the spherical system can be expressed as 



M6 v ) = ^ = - 3Sv{ \ +6) . (21) 

This expression is valid both in the linear and non-linear regimes and we refer to / as the total growth rate or non-linear 
growth rate. Along the Zeldovich curve, S v is a function of 8 and the growth rate can be expressed solely as a function of 8. 
Figure [5] shows the growth rate vs. (1 + 8) along the Zeldovich curve for different values of w at three different redshifts. The 
colour coding is same as that in figure [2] (1 + 8) is a more natural variable than 8 since the former is a measure of the total 
density of the system. As expected, the growth rate is higher for higher values of l + <5; larger densities grow faster, a signature 
of gravitational instability. At intermediate epochs, the growth rate is higher for lower values of w because the values of Q m 
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Figure 5. Non-linear growth factor as a function of (1 + 5) for models with four different values of to and the f! m = 1 model. Colour 
coding is same as in figure [2] Growth rates for the various models deviate only at intermediate epochs because the values of f! m are 
different in the different models. Lower value of w implies higher value of Q m and hence a higher growth rate. 



are higher (see figure [I]). The cumulative effect of higher growth factors can be seen by referring back to figure [3] Consider 
one of the collapsing perturbations, for example the one shown by the red dot in the upper and lower panels. It starts at the 
same point in the phase space at a — 0.1 for both models. At a = 1, it has evolved further along the Zeldovich curve for the 
fl m = 1 case as opposed to the ACDM case. The perturbations shown by the purple and violet dots are already out of the 
plot at a = 1 in the upper panel whereas are still at <5 < 5 in the lower panel. 

Using the fits in the above section, the non-linear growth rate in the regime — 1 ^ 5 ^ 1 is 



f 



3A[(1 + 8) - (1 + S) 1 



(22) 



Substituting for A and B using eqs. ( 15 1 and ( 16 1 and in the limit of small S, we get 



71+72 



1 + 



3 a L ra 



(23) 



To lowest order this reduces to the linear growth rate fu„ 
5. 



71+72 



At the next order, the non-linear correction is linear in 



5 CONNECTION TO OBSERVABLES: GALAXY CLUSTER PROFILES 



Studies of the matter density distribution and pattern of infall velocities around a rich cluster of galaxies have traditionally 
been used to provide constraints on the matter density Q m ( |Gunn fc Gottp972| 



Peebles||1976"l|Regos fc Geller||1989| . More 



recently, galaxy cluster profiles have also been used as a test of modified gravity ( LombTiser**eTa!*|2*oT2 1 . Here we investigate 
dependence of the cluster density and velocity profiles on the dark energy equation of state. We use the spherical infall model 
(e.g., |Silk fc Wirson||1979a||Silk fc Wilson||1979b| |Villumsen fc Davis] 1986||Lilje fc Lahav||1991| hereafter LL91) to track the 



evolution of the cluster into the mildly non-linear regime. Our overall set-up is similar to that described in Lahav et al. ( 1991 1 



5.1 Equations for spherical infall 

The initial configuration of the system is similar to that of the top-hat, but with a radial density variation and hence the 
perturbation cannot be described by a purely time dependent scale factor. The detailed description of the system and derivation 
of the relevant equations is given in Appendix [X] Here we state the main results. The system is described by a continuum of 
shells with Lagrangian coordinates X — r(ti)/ai, where r(ti) is the initial physical radius of the shell. The radius of the shell 
at any later time 

r(X,t) = b(X,t)X. (24) 
The initial perturbation is completely specified by its initial density and velocity profiles. Equations |2| and |3| generalise to 

5{X,U) = 1, 

pm (ti ) 

X (X t \ - - 1 ! 
5v[X ' U) ~ IfiHX^ti) ~ 

and the evolution of b(X, t) is given by 
"n ro ,i(l + A(X, U))a\ 



Hi 
2 



+ (l + 3to)fi w 



3(l+m) 



(25) 
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with initial conditions b(X,ti) = Ui, b(X,ti) = di(l + 8 v (X,ti)) and 

r-X 



A(X,U 



X 3 



6(X',U)X' 2 dX', 



where the spherically averaged density inside any general radius 'r' is given by 



AM = 3 



8{r')r' 2 dr' . 



The Jacobian that relates the initial Lagrangian and the physical coordinates is 



J{X,t) 



r dr 
X^dX 



1 + 



X db 



h a) 



and the density at any later time is 



S(X,t) 



{l + S(X,U))a 3 (1 + S(X, tj))a 3 



(26) 



(27) 



(28) 



(29) 



Note that the evolution of 5 is different from that given by eq. (JsJ) due to the spatial dependence of the perturbation scale 
factor b(X, t) and the combination 8 — 5 V no longer satisfies the Zeldovich relation. Instead the relevant quantity is the average 
density A. Substituting eq. (24 1 in eq. (271 and using the relations eqs. (26 1, (281 and (291 gives 

(1 + A{X,U))a 3 



A(X,t) 



b 3 



(30) 



Comparing eqs. ( |25[ ) and ( |30[ ) with eqs. Q and |8|, it is clear that the role played by 8 in the uniform density case is now 
played by A. Physically, this is expected because the velocity is set by the net mass inside the shell. The velocity divergence 
parameter is similar to eq. j9J, but with a spatial dependence 



S v (X,t) 



1 b(X,t) 



H b(X, t) 

The Zeldovich curve becomes a A — 8 V relation. 



(31) 



5.2 Initial conditions 

The initial density profile 8(X,ti) was chosen according to the prescription described in LL with some minor modifications. 



This prescription uses the theoretical framework of Bardeen et al. ( 1986 1 to predict the density profile around a primordial 



density peak that nucleates a cluster. The cosmological parameters today were set as Q m ,o = 0.29, f^.o = 0.71, Ho = lOO/i 
kms _1 Mpc _1 . The r.m.s. fluctuation on 8/i _1 Mpc, ag, fixed the normalisation of the initial power spectrum. We ran two sets 
of runs with erg = 0.8 and 0.9. In each set, three equation of state parameters were considered: w = —1,-1/2,-3/2. The 
linear density profile today was chosen to be the same for all three cases and the initial profiles at z = 1000 are obtained by 
multiplying by the appropriate linear growth rate factor. Therefore, the amplitude of the initial density profile was different 
for the three cases. The initial velocity profile was determined by imposing the linear Zeldovich condition in the A — 8 V phase 
space: 5 V (X,U) = -l/3fi/„ 1+72 A(X, U). 

Equation ( |25[ ) for b(X , t) was integrated numerically for a discrete set of shells from the initial time Z4 = 1000 to the 
final time Zf =0. During the course of evolution if b = 0, the shell is considered to have collapsed, and forms a part of the 
cluster core Q The initial set up is such that the inner shells collapse before the outer shells and no shell crossing takes place 
anywhere other than at the centre. Further details regarding the initial set-up can be found in Appendix [C] In the entire 
discussion, we assume that light traces matter exactly i.e., the bias factor between galaxies and dark matter is unity. The final 



density profile is computed using eq. ( 29 1 and the final infall velocity is 

v inf aii(X,t) = Hr(X,t) - r(X,t). (32) 
8 V and Vjn/aii are related as 

W) = _?^M (33) 
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Figure 6. The density and infall velocity profiles at z = for w = —1 (red, dotted) w = — 1/2 (green, dotdashed) and w = —3/2 
(brown, long dashed) cosmologies. The upper panel (lower) panels are w.r.t. X (r). Given the same linear amplitude today, the higher 
w models have higher amplitude in the non-linear regime (upper panel) because of differences in the growth factors. The w dependence 
gets weakened when plotted against the physical radius. 



5.3 Results 

Figure [6] shows the final evolved density and infall velocity profiles at z = for three cases: w = — 1 (red, dotted), w — —1/2 
(green, dotdashed) and w = —3/2 (brown, long dashed). The runs here correspond to initial conditions with erg = 0.8. The 
plots in the upper panel are vs. the Lagrangian coordinate (X) and those in the bottom panel are vs. the physical distance from 
the cluster centre (r). In the upper panel, the differences in evolution for different w models become apparent in the non- linear 
regime and the curves merge in the linear regime. This is because the linear theory power is fixed to be the same in all models 
at z = 0. Consequently, the amplitude of the initial profile at recombination is higher for higher w models because the growth 
factors are smaller and hence the final non-linear amplitude is more in these models. This is in qualitative agreement with 



results of McDonald, Trac, & Contaldi ( 2006 1, who examined the effect of w on the matter power spectrum. The sensitivity to 
w is almost eliminated when the curves are plotted against the physical radius r. In the absence of perturbations, the physical 
radius is just the comoving radius multiplied by the expansion factor r = a(t)X. For collapsing perturbations the net physical 
radius is always less than this; the higher the density and infall velocity, the more tightly bound the system and larger the 
difference. This causes a bunching of the curves when plotted against r. 

The situation is not too different at higher redshifts. At higher redshifts the distance from the cluster centre r (metric 
separation) is not the relevant observable and is replaced by the angular separation 9. For a flat universe, the two are related 
as 

(34) 



ax{a) ' 



x{ a ) is th e comoving distance to the object defined as 
da 



X (a) = c 



2 H(a)' 



(35) 



1 Ideally, for collisionless infall, a shell that collapses, bounces back and forth before eventually settling down to a periodic bouncing 
motion with maximum radius equal to 0.8 times its turnaround radius (jBertschinger 1985). In this case, the shell will collide with outer 
shells that are undergoing their initial collapse and the equations of motion will be differ from eq. l |-'5| . In this paper we follow [Lahav| 



et al. 



1 1991 ) and ignore these effects of secondary infall. 
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Figure 7. Density and infall velocity profiles at redshift z = 0.67. At higher redshifts, metric distance is replaced by the angular 
separation, which itself depends on the background cosmology. Plotting the profiles against the angular separation shows a increased 
deviation between models. However, isolating peculiar velocities at these redshifts is not possible with current observations. 



where H(a) is given by eq. |5| and c is the speed of light. Figure [7] shows the density and infall velocity patterns at z = 0.67 
plotted against the metric separation (top panel) and the angular separation (bottom panel). As is the case at z = 0, the 
density and velocity profiles when plotted against the metric separation, are not too sensitive to the change in w. But the 
angular separation also depends on w. For the same metric separation, higher w implies a larger angular distance. This 
broadens the curves and marginally improves the sensitivity to the equation of state parameter especially in the outer regions 
of the cluster. 

Determination of peculiar velocities requires having redshift independent distance measures. While such measures (e.g., 
Tully-Fisher relation) are well developed at low redshifts, their extension to higher redshifts is not well established. Isolating 
individual infall velocities at high redshift is extremely challenging and only statistical information is available currently. 
Ideally, the most relevant quantity for current observations would be the total radial velocity (r = Hr — v in f a u) and the 
corresponding zero- velocity surface (where r = 0). However, the zero- velocity surface was practically identical for the different 
values of w considered here (plot not shown), in agreement with the effect pointed out in Lahav et al. (19911. Figure [6] 
indicates that at z = and r ~ 3h~ Mpc, the deviation of the w = —1/2 curve from the w = —1 curves is about ~ 80 km/s 
(a deviation of about 13%) and figure [7] shows that at z = 0.67 and 9 ~ 10 arcmin the deviation is ~ 100 km/s (a deviation 
of about 25%). These deviations are just barely within the error bars for individual measurements ( Courtois fc Tuliy]|2012[ ) 
in the local universe, and probably beyond the current velocity resolution at high redshifts, but hopefully future advances in 
observational techniques will enable these measurements. 



5.4 Parameter degeneracies and the phase space picture 

The theoretically predicted profiles were obtained using fixed values of input parameters such as as (chosen to be 0.8). 
Determination of as has its own observational uncertainties which introduces a degeneracy between erg and w. Thus, using 
only the density or only the velocity information cannot not constrain w uniquely. We demonstrate that, at intermediate 
redshifts, this degeneracy can be removed by combining the density and velocity information by plotting it on the A — 5 V 
phase space plot. 

Figure [8] shows two sets of average density and velocity profiles and the phase space portrait at z = 0.67. The lower blue 
set corresponds to w = — 1 with as = 0.8 (dotted) and erg = 0.9 (dashed) and the upper red set is for w = —1/2 with the 
same erg values. The velocity profiles for w = — 1,(78 = 0.9 almost coincide with w — —1/2, as = 0.8. The same is true for the 
averaged density A. But, when the velocity information is combined with the density information and plotted on the phase 
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Figure 8. Phase space plot as a tool to lift the degeneracy of parameters. Blue (red) lines are for w = —1 (ui = —1/2) and dotted 
(dashed) lines are for erg = 0.8 (cr = 0.9). A higher value of <r has the same effect as a higher value of w and plotting only the infall 
velocity (first plot) or the spherically averaged overdensity (second plot) cannot constraint w. However, plotting the same information 
on the phase space portrait removes this degeneracy. The Zeldovich curve, which is a relation between the spherically averaged density 
A and fractional Hubble parameter 5 V , distinguishes cosmologies with different values of w. 



space portrait, the w = — 1 and w = —1/2 points clearly separate and practically lie on their corresponding Zeldovich curves. 
The relative deviation from the Zeldovich curve was calculated for all four sets of points and the maximum over all the points 
was 1.6 x 10~ 5 . 

The attracting nature of the Zeldovich curve ensures small uncertainties in the initial conditions do not affect the final 
state of the perturbations. Therefore, any degeneracy parameter that affects the initial conditions but not the evolution 
equations can be eliminated by plotting the density and velocity information on the phase space portrait. This can only work 
at intermediate redshifts (z ~ 1) where the Zeldovich curves are different for different w models. In addition, degeneracy 
parameters such as galaxy bias, which arise from observational constraints, cannot be eliminated by this plot. 



6 SUMMARY AND CONCLUSION 

We have obtained the non-linear density-velocity divergence relation by examining the dynamics of the spherical top-hat 
perturbations in the joint S — 5 V phase space. Imposing the 'equal age' condition on the perturbation parameters traces out a 
special curve in the instantaneous phase space of the system. Because of the non-autonomous nature of the S — S v evolution 
equations, this curve changes with time. We called this the Zeldovich curve and demonstrated that it acts like an 'attractor' 
for the dynamics of the perturbations. Small amplitude perturbations, as they evolve through the phase space asymptotically 
approach this curve and perturbations that start along the curve stay along it with a very high accuracy. Thus, the Zeldovich 
curve gives the long term behaviour of density and velocity perturbations and is the desired extension of the density- velocity 
curve in the non-linear regime. 

We obtained a 3% fit to the Zeldovich curve in the range — 1 S ^ 10 by generalising existing formulas of B92 and BC08. 
We found that, as in the linear regime, the non-linear relation mainly depends on the matter density parameter Q m and the 
explicit dependence on w is very weak. Although the result itself is not so surprising, the main new feature of this work is 
the interpretation of the relation in the phase space picture and the generality of the method. This work was restricted to 
flat cosmological models with dark matter and non-clustering dark energy described by a constant equation of state, but the 
approach described here can be extended to more complicated models of dark energy or alternate theories of gravity. 

Although the explicit dependence on w is weak, there is an implicit dependence through the evolution of Q m . As a practical 
application, we considered the evolution of density and velocity profiles of galaxies in the quasi- linear regime (8 ~ 5 — 10), 
before effects of virialization become important, and investigated the dependence on w. We found that the deviations in cluster 
profiles due to a 50% change in w are about ~ 100 km/s which are just marginally within the current observational errors 

. We also demonstrated that, at intermediate redshifts (z ~ 1), the density or 
velocity profiles alone cannot constrain w since these measurements are subject to degeneracies with parameters such as a%. 
However, these degeneracies can be eliminated by combining the information onto the A — 5 V phase space plot. This is because 
of the 'attracting' nature of the Zeldovich relation which makes it insensitive to small changes in the initial conditions. This, 
of course, does not eliminate the dependency on the galaxy bias parameter, which in this work was assumed to be unity. 

The Zeldovich curves were obtained using a spherical top-hat model. There are many sources of deviation from spherical 
symmetry in real systems; tidal interactions give rise to transverse forces, infall of smaller systems superimposes a random 
velocity on the radial infall ( |Diaferio fc G cllcr 1997) etc. Whether or not real data obey the Zeldovich relation even in 
a spherically averaged sense is not yet examined. Numerical simulations will be required to give more refined theoretical 
answers. On the observational side, isolating peculiar velocities at high redshift remains the main challenge since there are no 
well characterised redshift independent distance indicators at those redshifts. Nevertheless, we hope that the Zeldovich curve 



in the local universe (|Courtois fc Tully|[2012| 
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and its significance in the density-velocity phase space dynamics can be exploited to constrain cosmological parameters with 
future data. 
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APPENDIX A: EVOLUTION EQUATIONS FOR THE SPHERICAL PROFILE 



This appendix gives the equations that govern the evolution of a spherically symmetric matter perturbation in a background 
cosmology comprising of dark matter and dark energy described by a constant equation of state w. Dark matter and dark 
energy are uncoupled; dark matter clusters while dark energy stays spatially uniform. The equations derived here hold for 
perturbations with arbitrary radial density and velocity profiles. The spherical top-hat of [J2] becomes a special case. 

The physical configuration of the compensated spherical perturbation consists of an inner sphere surrounded by a com- 
pensating spherical shell. The background extends beyond the outer edge of the compensating shell. The density in the 
compensating shell is adjusted so that the mean density of the total mass enclosed within the perturbation is the same as 
that of the background. Let p m (p m ) denote the density of the background (perturbation). The origin is at the center of the 
inner sphere and r and x denote the physical and comoving distance from the center. 

We start with the generalised Euler equation for the evolution of the velocity ( Lima et al.|1997 Abramo et al.|2007 Pace 
et al.|2010| > 

du . „ ..<•' ^ . (A1) 



,, + u • V^u + 2-u 

at a 



= -a z S7 x <f>, 
= 4nGa 2 Sp m (^,t), 



(A2) 



where x is the comoving coordinate, u = dx/dt, <& is the peculiar gravitational potential and 5p m is the matter density 
perturbation. There are two main frameworks to describe the evolution of a fluid: Eulerian and Lagrangian. In the Eulerian 



description, one generally solves eq. (All for u as a function of some fixed grid coordinates, whereas in the Lagrangian 



description, the position is the main variable and is solved as a function of some initial coordinates and time. Since the main 
aim is to solve for the scale factor, the latter approach is more convenient. Combining the first two terms into a total derivative 



(d/dt + u • Vz = d/dt), using spherical symmetry and changing to physical radial coordinates r = ax, eqs. (All and (A2l 
become 



r r — — V r $, 

a 

where the dot is derivative w.r.t. time t and the background scale factor evolves as 



Hi 



j2 ti 



+ (i + 3uon*,i 



3(1+™) 



(A3) 
(A4) 

(A5) 



Define the Lagrangian coordinate of a shell as the initial comoving coordinate of the shell i.e., X — rj/oj. The physical 
radius of any shell at a later time can be written as 

r(X,t) = b(X,t)X, (A6) 

where b(X, i) can be thought of as the scale factor of the shell at X. By definition b(X, ti) = Oj. The perturbation is described 
by two quantities: initial overdensity parameter 

p m (X,U) 



S(X,ti) = 



(A7) 



and the velocity perturbation parameter 

-i^" 1 ' (A8) 

where f(X,ti) and p m {X,ti) are the initial total velocity and initial total density profiles of the perturbation and p m> i is the 
initial total density of the background. Conservation of mass implies that the matter density at any later time 

p m (X,ti)J(X,ti) 



' m{x ' t] J(X,t) 

where J is the Jacobian factor relating the Eulerian and Lagrangian volume elements, 
r 2 dr = J(X,t)X 2 dX. 



Using eq. ( A6 1, 



By construction, 
J(X,ti) = af. 



(A9) 
(A10) 
(All) 
(A12) 
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The background density p m — p m ,ia\lc? and using eqs. (A7l, |A9l and (A12l the perturbed density evolves as 

s Pm (x,t) - y — — 

From spherical symmetry, the force 



dr 



8pm(r )r dr 



Convert to Lagrangian coordinates by substituting eqs. (A6l, (AlOl and (A13l in (A14l 

„3 rX 



Define 



^Gp m ,iajbX 
3a 3 



3a A 
b 3 X 3 



(l + S(X',U))X' 2 dX\ 



X 3 



(8(X' y U))X ,2 dX'. 



(A13) 



(A14) 



(A15) 



(A16) 



A(X,ti) is the average fractional density of mass within radius X. Substituting eqs. (A6l and (A15l in eq. (A3 1, setting 
47rGp mil /3 = l/2H?n m ,i and using eq. (|A5| for the evolution of the background gives 



Hi 
2 



n m , t aUi + A(x,u)) 
b 3 



+ (i + 3w)n^ (^j 



3(1+™) 



(A17) 



The initial conditions are b(X, U) — at and b(X, ti) = di(l + S V (X, ti)). 

If the initial density is uniform, A{X,U) = Si and 8 V (X, U) = 8 v ,i and the X dependence drops out of the equation 
for evolution of b(X,t). An initially uniform spherical perturbation continues to stay uniform and can be described by 
r{X,t) = b(t)X. 



2 



b 3 



3(1+™) 



(A18) 



with initial conditions b(ti) = a% and b(U) = dj(l + 8 V) i). This equation is similar to many others stated more directly in the 



literature (see for e.g., Percival 20051. Here we use the route advocated in recent papers (Pace et al. 2010 Wintergerst & 



Pettorino|2010 1 so that the set up can be generalised to cases where dark energy may have more complicated behaviour. 

In the above derivation we have assumed that the Lagrangian system is a good coordinate system to describe the 
evolution. This is true as long as no shell crossing occurs within the system, which ensures that the mapping between the 
physical coordinates r and the Lagrangian coordinates X is unique. For the top-hat case, the argument that supports this 
assumption is presented in Appendix A of NC. Here we extend it to the spherical infall model. Let r e< j ge ,i be the initial 
physical distance of the edge of the cluster from the center. The velocity of the edge is r e d ge ,i = (1 + 5 v (X e d g e,ti))Hir e d ge ,i. 
Given the perturbation parameters S(X e d ge ,ti) and 8 v (X e d ge ,ti), one can always choose r e d ge ,i (and hence r e d ge ,i) arbitrarily 
small so that the time for the edge to reach any physical distance is arbitrarily large. The initial density 6(X, ti) is maximum 
at the center and decreases monotonically as X increases (see Appendix [c] for the initial set-up). The velocity perturbation 
S V (X, ti) is proportional to 8(X, ti), but with a negative sign and hence (1 + 8 V ) increases as one moves from the center to the 
edge. Thus, the net gravitational binding is tighter near the cluster center than the edge. This ensures that a shell initially 
closer to the center does not cross any shell at a greater initial radius and collapses into the center earlier. The choice of the 
edge radius fixes the net mass inside. The additional mass needed to satisfy mass conservation is put in a compensating shell 
between the cluster edge and the inner edge of the background and set on a critical trajectory outward so that it moves along 
with the background. 



APPENDIX B: TIME DEPENDENCE OF THE ZELDOVICH CONDITION FOR DARK ENERGY 

The evolution of the background and perturbation scale factors for the top-hat perturbation is given by 

£ _ _Hf fn m ,,a 3 _ ,„ /a,\ 3 ( 1+ra > 

a 2 



+ (1 + 3w)Sl*,i ( =^J J (Bl) 

/, - , \ /, ,(l + 3»)n*,(f) ). (B2) 
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where Qm,i, O^j are the Hubble parameter and density parameters at the initial time U and are related to their values today 
(a = 0) through 

ftm,i = _ 3 m '° - — _ 3(1+w y (B3) 

O „-3(i+ m ) 

n *>* - o =i(i+5o- ( B4 ) 

The initial conditions are a(ti) = a,, d(fi) = HiCii, b(ti) = a,i and b(ti) = HjOj(l + 5 Bl j). Substitute a; = a/a,i,y = &/<Xj, r = t^J/i 
in eqs. (jBlJ and ( |B2[ ). The two equations then read 

d^ = "2 ( v ^ + (1 + 3w) ^4xJ J' (B5) 

^ - ^ + (l + 3 W )n„y j, (B6) 

with initial conditions s(ri) = 1, i^Ti) = 1, y(ji) = 1, y(n) = (1 + &v,i), where the dots now denote the derivative w.r.t. 
t. Note that the equal age condition is unchanged by this scaling since the time is scaled by the same constant both for 
the background and perturbation. It is also clear that the equal age condition will only involve the parameters fi m ,i,fl</>,i 
and w. For a pure matter universe = and fl m ,i — 1 at all times. The ages of the background and perturbation are 



calculated analytically and the equal age condition reduces to eq. (Ill in the text. In the presence of a dark energy term this 



is not possible. For such cases, given a Si and <S„,j at r = r,, the age of the perturbation can be computed by integrating 



eq. (B6| back in time until y reaches zero. Fixing Si one performs a search in the S Vj i parameter space until the age of the 



background and perturbation are the same. However for models with w < — 1, an additional complication arises because 



eq. (B6l has a singularity when x = 0. This condition is independent of S v ,i and can complicate the search. So, instead of 
backward integration, we set y ~ at the big bang time (epoch when x = 0) and integrate forward in time until t;. The time 
elapsed gives the age of the perturbation. A search is performed to find the velocity at the bang time for which the ages of 
the background and perturbation are the same. The initial velocity at r, can be read off from the solution: S Vj t = y(Ti) — 1. 



APPENDIX C: SETTING THE INITIAL DENSITY PROFILE OF CLUSTER PROGENITORS 

The initial density and velocity profiles are set, as per the prescription in LL, with some minor modifications. This paper uses 
results from the BBKS ( Bardeen et al.|[l986 l analysis. The basic premise is that matter collapses around the peaks in the 
primordial density field, filtered on some appropriate length scale, forming the progenitors of bound cosmic structures. 
Start with the CDM power spectrum P(k) = PokT 2 (k), where the transfer function 

ln(l + 2.34q) 



T{q) 



[1 



■ + (16.1<?r + (5.46g) d + (6.71g) 4 



-1/4 



2.34<? 

q = fc/(fi m /i 2 Mpc _1 ). 

k is measured in /iMpc - . The normalisation Po is set by demanding that the r.m.s. fluctuation on 8 h' 1 Mpc (erg) defined as 



(CI) 
(C2) 



2 
°"8 



27T 2 



k 2 P(k) 



3(sin ka — ka cos ka) 
(kaf 



dk, 



where a = 8h 1 Mpc, matches the observed value. We ran separate cases with two values as = 0.; 
normalised power spectrum is filtered with a Gaussian filter of width 3/i _1 Mpc 

P f (k) = P(k)W 2 (R f ,k), 

where 

W(R f ,k) = e -^ R f k)2 . 



(C3) 

and erg = 0.9. The 
(C4) 

(C5) 



The height of any peak in the density distribution is usually denoted as vao where oq is the r.m.s. of the filtered density field. 
The ensemble average of the radial distribution of the density field around a peak of relative height v is given by 



S(r,u) = 



2tt 2 (to 



sin kr 
kr 



v — 7 2 !/ — 7# 



1 



T 



+ 



8Rl 



3 7 (l-7 2 ) 



dk. 



(06) 



Here 7 and _R* are spectral parameters and 9 is a function of v and 7 (see LL91). 

The fractional peak height v is a random variable with a comoving differential number density N p k{v) and ideally one 
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must evolve many initial profiles with different values of v drawn from this distribution. Instead, we follow LL91 and substitute 

e values (u) and 



for v and in eq. ( C6 1 the average values (v) and (6) 

r 



f°° vN vk (y)dv 
J N pk (y)dv 



J° °Af P k(y)dv 



(C7) 



where the threshold v± is the peak height above which clusters can form. This is set by requiring that the theoretical number 
density of galaxy clusters above the threshold equals the observed number density n p k,obs- 



n pk {vt) = / M p k{v)dv. 



(08) 



We use n 



pk.obs 



lQ~ b h 6 Mpc~ 6 for clusters with mass greater than 10 M e (Allen, Evrard, & Mantz 20111. 



The function <5(r, (u)) set up by this method corresponds to the linear profile at a ■ 
at any other redshift by multiplying by the appropriate growth factor. 

5(r, ( V )) a=ai = 5(r, (u)) a=1 D{a = ai) 



1 . This is transformed into the profile 



(C9) 



n , x 5 H(a) 
D(a) = -n m ,o^r^ 



(CIO) 



D(a= 1) ' 

where the growth factor D(a) is ( |Heath||1977| |Dodelson|2003 ) 

da' 

(a'H(a')/H ) 3 ' 

This expression for the growth factor has been derived for pressureless matter universes, but is also valid for dark energy 
cosmologies described by constant w. 1 5.1 sets up the initial profile in terms of Lagrangian coordinates: S(X, ti) = 8(r, {v)) a=ai ■ 
The initial profiles are chosen so that at z — the linear theory profile is same in all models. Since the growth factor is different 
for different w, the initial profiles have different amplitudes; higher w values have slower growth i.e., D(a — l)/D(a = ai) is 
smaller and hence a larger initial amplitude. The initial average density perturbation inside a radius X is 

(■X 



A(X,ti) = 



X 3 



X' 2 S(X' ,U)dX' . 



The velocity perturbation is chosen by imposing Zeldovich initial conditions 
1. 



« (x,t i ) = --n 



°- 5 ?A(X,U 



(Oil) 



(C12) 



At recombination, Q m ,i ~ 1 for all values of w. Note that the definition of 5 V slightly differs from a similar velocity perturbation 
parameter a defined in LL91. 
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